#ifndef LIKELIHOOD_HPP
#define LIKELIHOOD_HPP

#include "../Density.hpp"
#include <armadillo>
#include <iostream>
#include "../utils.hpp"
#include "Probitst.hpp"

using namespace arma;
using namespace std;

namespace Likelihoods{





class Probit : public Density::Tempering
{
		friend class Density::Tempering;
	public:
		Probit(double *Y, mat X) : Density::Tempering(Y, X){}
		~Probit();
		double Likelihood(mat theta)
		{
			double L=0;
			double sum=0;
			int m=theta.n_rows;
			mat bar=theta;//;.col(0);
			//Xb:
			for(int i=0;i<Get_n();i++)
			{	
			//	cout << theta.col(0) << "\\";
				mat foo= Get_X().row(i);
				sum=dot(foo,theta);
				
				if(Get_Y()[i]==1){
					L+=log(Phi(sum));
				}else if(Get_Y()[i]==0){
					L+=log(1-Phi(sum));
				}else{
					cout << "Warning values of Y!= {0,1}";
				}

			}
			//cout << exp(L);
			return L;
		}
		double Prior(mat theta){
			double res=0;
			return res;
		}

	
	/*private:
		double _phi_n;
		double _phi_n1;
		mat _X;
		double *_Y;
		int _p;
		int _b;
		int _n;*/
};

class Probit2 : public Density::IBIS
{

	friend class Density::IBIS;
	public:
	Probit2(double *Y, mat X, mat s) : Density::IBIS(Y,X){
			_s=s;
			_is=inv(s);
	}
	double	Density( mat theta, int n, int p)//pass the matrix as const
	{
		double L=0;
		double sum=0;
		int m=theta.n_rows;
		mat bar=theta;//.col(0);
		//Xb:
		for(int i=0;i<=p;i++)
		{	
		//	cout << theta.col(0) << "\\";
		//	cout << "test: " << n+i << "\n";
			mat foo= Get_X().row(n+i);
			sum=0;
			sum=dot(foo,theta);
			if(Get_Y()[n+i]==1){
				L+=log(Phi(sum));
			}else if(Get_Y()[n+i]==0){
				L+=log(1-Phi(sum));
			}else{
				cout << "Warning values of Y!= {0,1}";
			}
	
		}

	
		//cout << exp(L);
		return L;
	}

	double Prior(mat theta){
			
			double res=0;
			mat foo=theta.t();
			int d=theta.n_elem;
			mat bar=foo*_is*foo.t();
			//cout << "//" << bar;
			res=-d*0.5*(log(2*PI))-0.5*log(det(_s))-0.5*as_scalar(bar);

			return res;
		}
	inline mat Get_v(void){return _s;}

private:
	mat _s;
	mat _is;
	/*	mat _X;
	double *_Y;
	int _p;
	int _b;
	int _n;
*/
};


class Probit3 : public Density::GeomBridge
{
		friend class Density::GeomBridge;
	public:
		Probit3(double *Y, mat X, Distribution::Distribution *S,mat s) : Density::GeomBridge(Y, X, S){ 
			_s=s;
			_is=inv(s);
			d= new boost::math::normal_distribution<>(0,1); 
  }
		~Probit3()
		{
			delete d;
		}
		inline double Likelihood(mat theta1)
		{
			mat theta=theta1.t();
			double L=0;
			double sum=0;
			int m=theta.n_cols;
			//Xb
			int n=Get_n();
			mat X=Get_X();
			double *Y=Get_Y();
			for(int i=0;i<n;i++)
			{	
			//	cout << theta.col(0) << "\\";
				mat foo= X.row(i);
				sum=dot(foo,theta);
				double phi=Phi(sum);
				double phi2=1-Phi(sum);
				if(Y[i]==1){
					L+=log(phi);
				}else if(Y[i]==0){
					L+=log(phi2);
				}else{
					cout << "Warning values of Y!= {0,1}";
				}
			}
			//cout << exp(L);
			return L;

		}
		inline double Phi(double x)
		{

			double t=cdf((*d),x);
			if(t==1)
			{
				t=0.9999999;
			}
			if(t==0)
			{
				t=0.00000001;
			}
			return t;
			//return boost::math::cdf<double, boost::math::policies::policy<> >(boost::math::normal_distribution(0,1),x);
		}

		double Prior(mat theta){
			double res=0;
			mat foo=theta.t();
			int d=theta.n_elem;
			mat bar=foo*_is*foo.t();
			//cout << "//" << det(_s);
			double de=det(_s);
			if(de>100000){de=100000;}
			res=-d*0.5*(log(2*PI))-0.5*log(de)-0.5*as_scalar(bar);
			return res;
		}
		mat GradLik(mat theta)
		{
			mat X=Get_X();
			double *Y=Get_Y();
			int p=X.n_cols;
			mat sum(p,1);
			sum.fill(0);
			int n=X.n_rows;
			double xb=0;
			for(int i=0;i<n;i++)
			{
				xb=dot(X(i,span::all),theta);
				if(xb!=xb)
				{
					cout << "cb3";
				}
				double gpxb=Phi(xb);
				double pxb=pdf(*d,xb);
				double foo1=pxb*(1/gpxb);
				double foo2=pxb*(1/(1-gpxb));
                                //cout << "foo1 "<<foo1;
                                //cout << "foo2 "<<foo2;

				if(Y[i]==1)
				{	
					mat temp=foo1*X(i,span::all);
					sum=sum+temp.t();				
				}else{
								
					mat temp=foo2*X(i,span::all);
					sum=sum-temp.t();				
				}
			}
			sum+=-as_scalar(arma::sum(theta,0))/_s(0,0);
			
			return sum;

		}
		boost::math::normal_distribution<>* Get_d(void){return d;}
	        mat Get_s(void){return _s;}	
	private:
		mat _s;
		mat _is;
		boost::math::normal_distribution<> *d;
	/*private:
		double _phi_n;
		double _phi_n1;
		mat _X;
		double *_Y;
		int _p;
		int _b;
		int _n;*/
};

} // end namespace
#endif
